from general_neuro.conductance_neuron import Conductance_neuron
from general_neuro.alpha_synapse import Alpha_synapse
from scipy.integrate import odeint
from scipy import linspace,arange,linalg,matrix,vectorize,array,hstack,vstack
from pylab import plot,show,subplot
import pylab

#----------------------------------------------------------------------------------
# some synaptic input
#----------------------------------------------------------------------------------
holding_V = [-79.0,-70.0,-50.0] # in mV
index = 0

# ----------------- simulate and get results for the first hp ------------
c = Conductance_neuron(g_leak = 3.0, E_leak = -56.5, tau = 100, V_init = -56.5)

c.register_input("E",Alpha_synapse(onset=14.5,tau=3.0,gmax=0.006,E_rev=10.0))
c.register_input("I",Alpha_synapse(onset=11,tau=8.5,gmax=0.0175,E_rev=-79.0))

t = arange(0,200,0.1)

current = []
for index in range(3):
    c.set_V_clamp(holding_V[index]) # set the neuron into V_clamp mode.

    sv_init = c.get_sv_init_list()
    sv_results = odeint(c.dX_dt, sv_init, t)

    # reconstruct the various currents from the sv_results
    current.append( [] )
    for i in range(len(t)):
        time = t[i]
        current[index].append( c.calculate_input_current("V_clamp",sv_results[i],time) )
    current[index] = array(current[index])

def solve_for_conductances(c1, c2, E_rev1, E_rev2, hp1, hp2):
    """
    c1 = g1 * (E_rev1 - hp1) + g2 * (E_rev2 - hp1)
    c2 = g1 * (E_rev1 - hp2) + g2 * (E_rev2 - hp2)

         |     c1         (E_rev2-hp1)|
         |                            |
         |     c2         (E_rev2-hp2)|
    g1 = ------------------------------
         |(E_rev1-hp1)    (E_rev2-hp1)|
         |                            |
         |(E_rev1-hp2)    (E_rev2-hp2)|

         |(E_rev1-hp1)          c1    |
         |                            |
         |(E_rev1-hp2)          c2    |
    g2 = ------------------------------
         |(E_rev1-hp1)    (E_rev2-hp1)|
         |                            |
         |(E_rev1-hp2)    (E_rev2-hp2)|
    """
    num_g1  = [[c1, E_rev2-hp1],[c2, E_rev2-hp2]]
    num_g2  = [[E_rev1-hp1, c1],[E_rev1-hp2, c2]]
    denom   = [[E_rev1-hp1, E_rev2-hp1],[E_rev1-hp2, E_rev2-hp2]]
    num_g1 = linalg.det(matrix(num_g1))
    num_g2 = linalg.det(matrix(num_g2))
    denom = linalg.det(matrix(denom))
    g1 = num_g1/denom
    g2 = num_g2/denom
    return (g1,g2)

def sfc(c1,c2,E_rev1,E_rev2,hp1,hp2):
    a = matrix([[E_rev1-hp1, E_rev2-hp1],[E_rev1-hp2, E_rev2-hp2]])
    b = matrix([[c1],[c2]])
    x = linalg.solve(a,b)
    return tuple(x)

def c3(c1,c2,v1,v2,v3,ve,vi):
    # from mathematica
    n1 = -(c2 * v1 - c1 * v2 + c1 * ve - c2 * ve) * (-v3 + vi)
    d1 = (v1-v2)*(ve-vi)
    n2 = -(c2 * v1 - c1 * v2 + c1 * vi - c2 * vi) * (-v3 + ve)
    d2 = (-v1+v2)*(ve-vi)
    return n1/d1+n2/d2
    
'''
# let's calculate the conductances from the currents.
g1 = []
g2 = []
print I[0].keys()
for i in range(len(t)):
    g1element,g2element = solve_for_conductances(I[0]["V_clamp"][i],I[1]["V_clamp"][i],\
                            10,-79,holding_V[0],holding_V[1])
    g1.append(g1element)
    g2.append(g2element)
'''

sfc_vectorized = vectorize(sfc)
g11,g21 = sfc_vectorized(current[0],current[1],-20.0,-79.0,holding_V[0],holding_V[1])
c3_vectorized = vectorize(c3)
c3_result = c3_vectorized(current[0],current[1],holding_V[0],holding_V[1],holding_V[2],-500,-79)
def c3_small(c1,c2,v1,v2,v3):
    n1 = c2 * v1 - c1 * v2 + c1 * v3 - c2 * v3
    d1 = v1 - v2
    return n1/d1

def ve(c3,v3,ge,gi,vi):
    if ge > 0.8*maxge:
        return (c3 + ge*v3 + gi * v3 - gi*vi)/ge
    else:
        return 0

ve_vectorized = vectorize(ve)
maxge = g11.max()
subplot(311)
pylab.title("From Currents to Conductances.")
plot(t,current[2])
pylab.ylabel("V_clamp at %2.2f"%(holding_V[2]))
subplot(312)
plot(t,g11,t,g21)
pylab.ylabel("conductances")
subplot(313)
plot(t,ve_vectorized(current[2],holding_V[2],g11,g21,-79.00000))
pylab.ylabel("ve")
show()

